0. パッケージの読み込み
本稿はRの基本操作とtidyverseパッケージによるデータハンドリングができることを前提としている。tidyverseパッケージを用いたデータ処理については、以下の書籍などを参照。
R for Data Science (Wickham & Grolemund, 2016)
電子書籍, 日本語R Graphics Coocbook 2nd Edition (Chang, 2018)
電子書籍, 日本語RユーザのためのRstudio[実践]入門~tidyverseによるモダンな分析フローの世界 改訂2版 (松村 et al., 2021) 出版社サイト
使用するパッケージは以下のとおりである。
## データハンドリング
library(tidyverse)
library(readxl)
library(knitr)
library(easystats)
library(DT)
library(scales)
library(emmeans)
## モデリング
library(brms)
library(rethinking)
library(rstan)
library(rstanarm)
library(cmdstanr)
library(ggeffects)
library(DHARMa)
library(DHARMa.helpers)
## ネットワーク分析
library(sna)
library(ANTs)
library(asnipe)
library(aninet)
library(vegan)
library(igraph)
library(tidygraph)
library(ggraph)
## グラフや表関連
library(ggtext)
library(ggforce)
library(ggbeeswarm)
library(plotly)
library(bayesplot)
library(viridis)
library(ggnewscale)
library(GGally)
library(flextable)
library(ggrepel)
library(patchwork)
library(kableExtra)
library(ggsci)
library(lemon)
library(ggsignif)
library(ggh4x)
library(ggrepel)
library(ggpattern)
# ## フォント関連
library(extrafont)
require(systemfonts)
# ベイズモデルの設定
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) 1 データの読み込み
以下、基本的なデータの読み込みを行う。なお、本論文で使用するのはた8つの期間（2018交尾期、2019非交尾期、2019交尾期、2020非交尾期、2020交尾期、2021非交尾期、2021交尾期、2022非交尾期）のデータである。
1.1 Daily data の読み込み
以下、日ベースのデータについてデータフレームを作成する。
なお、**には観察年の下2桁と交尾期(m)/非交尾期(nm)が入る。
base_**: 追跡時間や気温についての情報。
female_**: 各メスの確認の有無と発情の有無、年齢など（分派集団で確認された個体も含む）
male_**: 各オスの確認の有無
female_all: 全期間をまとめたメスの確認状況
male_all: 全期間をまとめたオスの確認状況
1.2 個体追跡データの読み込み
以下では、個体追跡の生データとして、以下のデータシートを読み込む。
focal_18m: 2018年交尾期の個体追跡データ
focal_19m: 2019年交尾期の個体追跡データ
focal_20m: 2020年交尾期の個体追跡データ
focal_21m: 2021年交尾期の個体追跡データ
focal_19nm: 2019年非交尾期の個体追跡データ
focal_21nm: 2021年非交尾期の個体追跡データ
focal_22nm: 2022年非交尾期の個体追跡データ
focal_combined: 全個体追跡データを結合したもの
2 メス同士のCSIの算出
2.1 CSIの定義
ここでは、TYとのCSIを算出する。
CSIの算出式は以下の通り。なお、\(i = 1,2,3,...,N\)はメスのIDを表す。\(P_i\)の算出の際には互いに毛づくろいしているポイントは除いた。対象としたのは2018年時点で6歳以上だった非発情の個体である。
- \(G_{ij}\): メスiとメスjの毛づくろい時間割合 (iとjが毛づくろいしたポイント数/iとj総ポイント数)
- \(P_{i}\): メスiとメスjの近接時間割合 (3m以内にいとjがいたポイント数/iまたはjが地上採食または休息したポイント数)
- \(\overline{G}: \frac{1}{N} \sum^{N}_{i=1}G_{i}\)
- \(\overline{P}: \frac{1}{N} \sum^{N}_{i=1}P_{i}\)
\[ CSI_{i} = \frac{\left(\frac{G_{i}}{\overline{G}} +\frac{P_{i}}{\overline{P}}\right)}{2} \;\;(i = 1,2,\dots,N) \]
使用するのは2018年交尾期、2019年交尾期、2019年非交尾期、2021年非交尾期、2022年非交尾期のデータである。
2.2 データの加工
交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_b、focal_19m_b)。
focal_18m %>%
left_join(female_18m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_18m_b
focal_19m %>%
left_join(female_19m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_19m_b2.3 CSIの算出(全期間)
算出するメスのIDは、全期間のデータがあった以下の通り。
femaleID <- sort(c("Kil","Mik","Koh","Aka","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))2.3.1 毛づくろい時間割合の算出
focal_all <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm, focal_21nm, focal_22nm)
## 毛づくろい相手を表す列を作成
focal_all %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
focal_all_b %>%
filter(subject %in% femaleID) %>%
group_by(subject) %>%
summarise(dur = n()) -> focal_duration対象となるメスと毛づくろいしていたポイントを抽出する。
focal_all_b %>%
filter(activity == "G") %>%
filter(subject %in% femaleID) %>%
select(date, no_focal, time, subject, groom, groom2) %>%
pivot_longer(cols = c("groom","groom2"),
names_to = "groom", values_to = "ID") %>%
filter(ID %in% femaleID) -> groom_pairs最後に、毛づくろい時間割合を算出してマトリックスにする。
2.3.2 近接時間割合の算出
続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。
まず、地上採食・休息・毛づくろいのポイントのみを抽出する。
focal_all_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>%
replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。
prox_denom <- matrix(nrow = length(femaleID), ncol = length(femaleID))
for(i in seq_along(femaleID)){
for(j in seq_along(femaleID)){
prox_denom[i,j] <- focal_prox %>%
filter(subject == femaleID[i]) %>%
filter(groom != femaleID[j] & groom2 != femaleID[j]) %>%
nrow()
}
}
## 転置行列と足す
prox_denom2 <- prox_denom + t(prox_denom)
## 対角成分は0にする
diag(prox_denom2) <- 0
rownames(prox_denom2) <- femaleID
colnames(prox_denom2) <- femaleID続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。
prox_numer <- matrix(nrow = length(femaleID), ncol = length(femaleID))
for(i in seq_along(femaleID)){
for(j in seq_along(femaleID)){
prox_numer[i,j] <- focal_prox %>%
filter(subject == femaleID[i]) %>%
filter(groom != femaleID[j] & groom2 != femaleID[j]) %>%
filter(str_detect(x0_1m, femaleID[j])|str_detect(x1_3m, femaleID[j])) %>%
nrow()
}
}
## 転置行列と足す
prox_numer2 <- prox_numer + t(prox_numer)
## 対角成分は0にする
rownames(prox_numer2) <- femaleID
colnames(prox_numer2) <- femaleID最後に、これらを基に近接時間割合を算出する。
2.3.3 CSIの算出
最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。
library(aninet)
CSI_female <- dyadic_csi(list(groom_mat, prox_mat))
colnames(CSI_female) <- femaleID
rownames(CSI_female) <- femaleIDCSIを図示すると以下のようになる。
CSI_female %>%
data.frame() %>%
rownames_to_column(var = "ID1") %>%
pivot_longer(2:16,
names_to = "ID2",
values_to = "CSI") %>%
ggplot(aes(x = ID1, y = ID2))+
geom_tile(aes(fill = CSI))+
scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
midpoint = 1)+
theme(aspect.ratio = 1)+
labs(x = "", y = "")
2.4 CSIの算出(2018交尾期～2019交尾期)
算出するメスのIDは、この期間のデータがあった以下の通り。
femaleID2 <- sort(c("Kil","Mik","Koh","Aka","Tam","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))2.4.1 毛づくろい時間割合の算出
focal_fh <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm)
## 毛づくろい相手を表す列を作成
focal_fh %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
focal_fh_b %>%
filter(subject %in% femaleID2) %>%
group_by(subject) %>%
summarise(dur = n()) -> focal_duration_fh対象となるメスと毛づくろいしていたポイントを抽出する。
focal_fh_b %>%
filter(activity == "G") %>%
filter(subject %in% femaleID2) %>%
select(date, no_focal, time, subject, groom, groom2) %>%
pivot_longer(cols = c("groom","groom2"),
names_to = "groom", values_to = "ID") %>%
filter(ID %in% femaleID2) -> groom_pairs_fh最後に、毛づくろい時間割合を算出してマトリックスにする。
2.4.2 近接時間割合の算出
続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。
まず、地上採食・休息・毛づくろいのポイントのみを抽出する。
focal_fh_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>%
replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox_fh続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。
prox_denom_fh <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))
for(i in seq_along(femaleID2)){
for(j in seq_along(femaleID2)){
prox_denom_fh[i,j] <- focal_prox_fh %>%
filter(subject == femaleID2[i]) %>%
filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>%
nrow()
}
}
## 転置行列と足す
prox_denom_fh2 <- prox_denom_fh + t(prox_denom_fh)
## 対角成分は0にする
diag(prox_denom_fh2) <- 0
rownames(prox_denom_fh2) <- femaleID2
colnames(prox_denom_fh2) <- femaleID2続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。
prox_numer_fh <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))
for(i in seq_along(femaleID2)){
for(j in seq_along(femaleID2)){
prox_numer_fh[i,j] <- focal_prox_fh %>%
filter(subject == femaleID2[i]) %>%
filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>%
filter(str_detect(x0_1m, femaleID2[j])|str_detect(x1_3m, femaleID2[j])) %>%
nrow()
}
}
## 転置行列と足す
prox_numer_fh2 <- prox_numer_fh + t(prox_numer_fh)
## 対角成分は0にする
rownames(prox_numer_fh2) <- femaleID2
colnames(prox_numer_fh2) <- femaleID2最後に、これらを基に近接時間割合を算出する。
2.4.3 CSIの算出
最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。
library(aninet)
CSI_fh_female <- dyadic_csi(list(groom_mat_fh, prox_mat_fh))
colnames(CSI_fh_female) <- femaleID2
rownames(CSI_fh_female) <- femaleID2CSIを図示すると以下のようになる。
CSI_fh_female %>%
data.frame() %>%
rownames_to_column(var = "ID1") %>%
pivot_longer(2:16,
names_to = "ID2",
values_to = "CSI") %>%
ggplot(aes(x = ID1, y = ID2))+
geom_tile(aes(fill = CSI))+
scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
midpoint = 1)+
theme(aspect.ratio = 1)+
labs(x = "", y = "")
2.5 CSIの算出(Tam生存前まで)
算出するメスのIDは、この期間のデータがあった以下の通り。
femaleID2 <- sort(c("Kil","Mik","Koh","Aka","Tam","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))2.5.1 毛づくろい時間割合の算出
focal_withTam <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm, focal_20m, focal_21nm) %>%
mutate(date = as_date(date)) %>%
filter(date <= "2021-03-14")
## 毛づくろい相手を表す列を作成
focal_withTam %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_withTam_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
focal_withTam_b %>%
filter(subject %in% femaleID2) %>%
group_by(subject) %>%
summarise(dur = n()) -> focal_duration_withTam対象となるメスと毛づくろいしていたポイントを抽出する。
focal_withTam_b %>%
filter(activity == "G") %>%
filter(subject %in% femaleID2) %>%
select(date, no_focal, time, subject, groom, groom2) %>%
pivot_longer(cols = c("groom","groom2"),
names_to = "groom", values_to = "ID") %>%
filter(ID %in% femaleID2) -> groom_pairs_withTam最後に、毛づくろい時間割合を算出してマトリックスにする。
2.5.2 近接時間割合の算出
続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。
まず、地上採食・休息・毛づくろいのポイントのみを抽出する。
focal_withTam_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>%
replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox_withTam続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。
prox_denom_withTam <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))
for(i in seq_along(femaleID2)){
for(j in seq_along(femaleID2)){
prox_denom_withTam[i,j] <- focal_prox_withTam %>%
filter(subject == femaleID2[i]) %>%
filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>%
nrow()
}
}
## 転置行列と足す
prox_denom_withTam2 <- prox_denom_withTam + t(prox_denom_withTam)
## 対角成分は0にする
diag(prox_denom_withTam2) <- 0
rownames(prox_denom_withTam2) <- femaleID2
colnames(prox_denom_withTam2) <- femaleID2続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。
prox_numer_withTam <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))
for(i in seq_along(femaleID2)){
for(j in seq_along(femaleID2)){
prox_numer_withTam[i,j] <- focal_prox_withTam %>%
filter(subject == femaleID2[i]) %>%
filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>%
filter(str_detect(x0_1m, femaleID2[j])|str_detect(x1_3m, femaleID2[j])) %>%
nrow()
}
}
## 転置行列と足す
prox_numer_withTam2 <- prox_numer_withTam + t(prox_numer_withTam)
## 対角成分は0にする
rownames(prox_numer_withTam2) <- femaleID2
colnames(prox_numer_withTam2) <- femaleID2最後に、これらを基に近接時間割合を算出する。
2.5.3 CSIの算出
最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。
library(aninet)
CSI_withTam_female <- dyadic_csi(list(groom_mat_withTam, prox_mat_withTam))
colnames(CSI_withTam_female) <- femaleID2
rownames(CSI_withTam_female) <- femaleID2CSIを図示すると以下のようになる。
CSI_withTam_female %>%
data.frame() %>%
rownames_to_column(var = "ID1") %>%
pivot_longer(2:16,
names_to = "ID2",
values_to = "CSI") %>%
ggplot(aes(x = ID1, y = ID2))+
geom_tile(aes(fill = CSI))+
scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
midpoint = 1)+
theme(aspect.ratio = 1)+
labs(x = "", y = "")
2.6 CSIの算出(調査期間ごと)
2.6.1 毛づくろい時間割合の算出
各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
focal_all_b %>%
filter(subject %in% femaleID) %>%
group_by(subject, study_period) %>%
summarise(dur = n()) -> focal_duration_sp対象となるメスと毛づくろいしていたポイントを抽出する。
focal_all_b %>%
filter(activity == "G") %>%
filter(subject %in% femaleID) %>%
select(date, no_focal, time, subject, groom, groom2, study_period) %>%
pivot_longer(cols = c("groom","groom2"),
names_to = "groom", values_to = "ID") %>%
filter(ID %in% femaleID) -> groom_pairs最後に、毛づくろい時間割合を算出してマトリックスにする。
groom_mat_sp <- list()
sp <- unique(groom_pairs$study_period)
for(i in c(1,2,3,5)){
groom_mat_sp[[i]] <- df.to.mat(groom_pairs %>% filter(study_period == sp[i]),
actor = "subject",
receiver = "ID",
sym = T,
tob = focal_duration_sp %>% filter(study_period == sp[i]) %>%.$dur)
}
## 2021非交尾期はホタルがfemaleIDに含まれるメスと毛づくろいをしなかったので、別に求める。
groom_mat_sp[[4]] <- df.to.mat(groom_pairs %>% filter(study_period == sp[4]),
actor = "subject",
receiver = "ID",
sym = T,
tob = focal_duration_sp %>%
filter(study_period == sp[4]) %>%
filter(subject != "Hot") %>%
.$dur)
## ホタルの行列を追加
groom_mat_sp[[4]] %>%
data.frame() %>%
bind_rows(data.frame(row.names = "Hot")) %>%
mutate(Hot = 0L)-> groom_mat_sp4
## 行名と列名を並べ替え
groom_mat_sp4[15,1:15] <- 0
groom_mat_sp[[4]] <- groom_mat_sp4 %>%
as.matrix()
groom_mat_sp[[4]] <- groom_mat_sp[[4]][sort(rownames(groom_mat_sp[[4]])),sort(colnames(groom_mat_sp[[4]]))] 2.6.2 近接時間割合の算出
続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。
続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。
mat <- matrix(nrow = length(femaleID), ncol = length(femaleID))
prox_denom_sp <- list(mat,mat,mat,mat,mat)
prox_denom_sp2 <- list(mat,mat,mat,mat,mat)
for(i in seq_along(sp)){
for(j in seq_along(femaleID)){
for(k in seq_along(femaleID))
prox_denom_sp[[i]][j,k] <- focal_prox %>%
filter(study_period == sp[i]) %>%
filter(subject == femaleID[j]) %>%
filter(groom != femaleID[k] & groom2 != femaleID[k]) %>%
nrow()
diag(prox_denom_sp[[i]]) <- 0
rownames(prox_denom_sp[[i]]) <- femaleID
colnames(prox_denom_sp[[i]]) <- femaleID
prox_denom_sp2[[i]] <- prox_denom_sp[[i]] + t(prox_denom_sp[[i]])
}
}続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。
prox_numer_sp <- list(mat,mat,mat,mat,mat)
prox_numer_sp2 <- list(mat,mat,mat,mat,mat)
for(i in seq_along(sp)){
for(j in seq_along(femaleID)){
for(k in seq_along(femaleID))
prox_numer_sp[[i]][j,k] <- focal_prox %>%
filter(study_period == sp[i]) %>%
filter(subject == femaleID[j]) %>%
filter(groom != femaleID[k] & groom2 != femaleID[k]) %>%
filter(str_detect(x0_1m, femaleID[k])|str_detect(x1_3m, femaleID[k])) %>%
nrow()
diag(prox_numer_sp[[i]]) <- 0
rownames(prox_numer_sp[[i]]) <- femaleID
colnames(prox_numer_sp[[i]]) <- femaleID
prox_numer_sp2[[i]] <- prox_numer_sp[[i]] + t(prox_numer_sp[[i]])
}
}最後に、これらを基に近接時間割合を算出する。
2.6.3 CSIの算出
最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。
library(aninet)
CSI_sp_female <- list(mat,mat,mat,mat,mat)
for(i in seq_along(sp)){
CSI_sp_female[[i]] <- dyadic_csi(list(groom_mat_sp[[i]], prox_mat_sp[[i]]))
colnames(CSI_sp_female[[i]]) <- femaleID
rownames(CSI_sp_female[[i]]) <- femaleID
}CSIを図示すると以下のようになる。
CSI_sp_df <- data.frame()
for(i in seq_along(sp)){
CSI_df <- CSI_sp_female[[i]] %>%
data.frame() %>%
rownames_to_column(var = "ID1") %>%
pivot_longer(2:16,
names_to = "ID2",
values_to = "CSI") %>%
mutate(study_period = sp[i])
CSI_sp_df <- bind_rows(CSI_sp_df, CSI_df)
}
CSI_sp_df %>%
ggplot(aes(x = ID1, y = ID2))+
geom_tile(aes(fill = CSI))+
scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
midpoint = 1)+
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = -45,
size = 7),
axis.text.y = element_text(angle = -45,
size = 7))+
labs(x = "", y = "")+
facet_rep_wrap(~study_period, repeat.tick.labels = TRUE)
各調査期間のCSIマトリックス間の相関は以下の通り。
3 TY・ITとメスのCSIの計算
3.1 CSIの定義
ここでは、TY・ITとのCSIを算出する。
算出方法は Yamaguchi & Kazahari (2022) に従う。CSIの算出式は以下の通り。なお、\(i = 1,2,3,...,N\)はメスのIDを表す。\(P_i\)の算出の際にはTY/ITと毛づくろいしているポイントは除いた。対象としたのは2018年時点で5歳以上だった非発情の個体である。
- \(G_{i}\): メスiとTY/ITの毛づくろい時間割合 (TYと毛づくろいしたポイント数/総ポイント数)
- \(P_{i}\): メスiとTY/ITの近接時間割合 (3m以内にTYがいたポイント数/地上採食または休息したポイント数)
- \(\overline{G}: \frac{1}{N} \sum^{N}_{i=1}G_{i}\)
- \(\overline{P}: \frac{1}{N} \sum^{N}_{i=1}P_{i}\)
\[ CSI_{i} = \frac{\left(\frac{G_{i}}{\overline{G}} +\frac{P_{i}}{\overline{P}}\right)}{2} \;\;(i = 1,2,\dots,N) \]
使用するのは2018年交尾期、2019年交尾期、2019年非交尾期、2021年非交尾期、2022年非交尾期のデータである。
3.2 TYとメスのCSIの算出
3.2.1 データの読み込みと加工
まずは、TYの出入りがあった2019年交尾期、2021年非交尾期についてTYが確認できた時間帯のデータを読み込む。
TY_19m <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
sheet = "male_presence_long") %>%
filter(maleID == "TY") %>%
mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
hour(.),minute(.))) %>%
mutate(date = as_date(date)) %>%
## 分派集団にいるときは除く
filter(groupID != 41 & groupID != 53 & groupID != 54)
TY_21nm <- read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
filter(maleID == "TY") %>%
mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
hour(.),minute(.))) %>%
mutate(date = as_date(date)) 上記2つの期間の個体追跡データに、各個体追跡セッション中にTYが確認できていたかを表す列を追加し、TYが確認できたものだけを抽出する(focal_19m_b、focal_21nm_b)。
focal_19m %>%
distinct(no_focal, date, start_time, fin_time) %>%
left_join(TY_19m %>% select(date, male_presence, first, last), by = "date") %>%
## 個体追跡セッション中にTYが確認されていたか
mutate(TY = ifelse((first <= start_time & last >= fin_time)|
(first >= start_time & first <= fin_time)|
(last >= start_time & last <= fin_time), 1, 0)) %>%
replace_na(list(TY = 0)) -> focal_19m_TY
focal_19m %>%
left_join(focal_19m_TY %>% select(no_focal, TY), by = "no_focal") %>%
filter(TY == "1") -> focal_19m_TY_b
focal_21nm %>%
distinct(no_focal, date, start_time, fin_time) %>%
left_join(TY_21nm %>% select(date, male_presence, first, last), by = "date") %>%
## 個体追跡セッション中にTYが確認されていたか
mutate(TY = ifelse((first <= start_time & last >= fin_time)|
(first >= start_time & first <= fin_time)|
(last >= start_time & last <= fin_time), 1, 0)) %>%
replace_na(list(TY = 0)) -> focal_21nm_TY
focal_21nm %>%
left_join(focal_21nm_TY %>% select(no_focal, TY), by = "no_focal") %>%
filter(TY == "1") -> focal_21nm_TY_b交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_TY_b、focal_19m_TY_b)。
focal_18m %>%
left_join(female_18m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_18m_TY_b
focal_19m_TY_b %>%
left_join(female_19m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_19m_TY_b3.2.2 CSIの算出(全調査期間)
3.2.2.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 全期間を結合
focal_all_TY <- bind_rows(focal_18m_TY_b, focal_19m_TY_b, focal_19nm,
focal_21nm_TY_b, focal_22nm)
## 毛づくろい相手を表す列を追加
focal_all_TY %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_TY_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のTYとの毛づくろい時間を算出する。
focal_all_TY_b %>%
filter(groom == "TY"|groom2 == "TY") %>%
group_by(subject) %>%
summarise(no_groom = n()) %>%
ungroup()-> no_groom_TY最後に、毛づくろい時間割合を算出する。
3.2.2.2 近接時間割合
地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。
focal_all_TY_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_TY分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。
TYと3m以内に近接していたポイント数を算出する。
focal_prox_TY %>%
filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>%
group_by(subject) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_TY最後に、近接時間割合を算出する。
3.2.2.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_TY <- mean(prop_groom_TY$prop_groom)
## 平均近接時間割合
mean_prox_TY <- mean(prop_prox_TY$prop_prox)
## CSIの算出
### 平均値
prop_groom_TY %>%
left_join(prop_prox_TY, by = "subject") %>%
mutate(CSI_TY = ((prop_groom/mean_groom_TY) + (prop_prox/mean_prox_TY))/2) -> CSI_TY算出した値が以下の通り。Kitは一日TYに長時間毛づくろいした日があり、それが強く影響したようだ。
CSI_TY %>%
mutate(subject = fct_reorder(subject, CSI_TY)) %>%
ggplot(aes(x = subject, y = CSI_TY))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(paste(CSI, " with ", italic(TY))))+
labs(title = "A")+
coord_flip(ylim = c(0,5.3)) +
scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
expand= c(0,0))+
theme_bw()+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "bold.italic",
hjust=0,
family = "Times New Roman"),
axis.text.x = element_text(family = "Times New Roman",
face = "bold"),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8,
family = "Times New Roman"),
axis.title.x = element_text(size = 10.5,
family = "Times New Roman"),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(family = "Times New Roman",
hjust = 0)) -> p_CSI_TY
p_CSI_TY
3.2.3 CSIの算出(2018交尾期～2019交尾期)
3.2.3.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 期間を結合
focal_fh_TY <- bind_rows(focal_18m_TY_b, focal_19m_TY_b, focal_19nm)
## 毛づくろい相手を表す列を追加
focal_fh_TY %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_TY_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のTYとの毛づくろい時間を算出する。
focal_fh_TY_b %>%
filter(groom == "TY"|groom2 == "TY") %>%
group_by(subject) %>%
summarise(no_groom = n()) %>%
ungroup() -> no_groom_fh_TY最後に、毛づくろい時間割合を算出する。
3.2.3.2 近接時間割合
地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。
focal_fh_TY_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_fh_TY分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。
focal_prox_fh_TY %>%
group_by(subject) %>%
summarise(dur = n()) %>%
ungroup() -> prox_duration_fh_TYTYと3m以内に近接していたポイント数を算出する。
focal_prox_fh_TY %>%
filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>%
group_by(subject) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_fh_TY最後に、近接時間割合を算出する。
3.2.3.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_fh_TY <- mean(prop_groom_fh_TY$prop_groom)
## 平均近接時間割合
mean_prox_fh_TY <- mean(prop_prox_fh_TY$prop_prox)
## CSIの算出
prop_groom_fh_TY %>%
left_join(prop_prox_fh_TY, by = "subject") %>%
mutate(CSI_TY = ((prop_groom/mean_groom_fh_TY) + (prop_prox/mean_prox_fh_TY))/2) -> CSI_fh_TY算出した値が以下の通り。
CSI_fh_TY %>%
mutate(subject = fct_reorder(subject, CSI_TY)) %>%
ggplot(aes(x = subject, y = CSI_TY))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(CSI[i]))+
coord_flip(ylim = c(0,5.3)) +
scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
expand= c(0,0))+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 10.5),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(hjust = 0.5)) 
3.2.4 CSIの算出(非交尾期のみ)
3.2.4.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 全期間を結合
focal_nm_TY <- bind_rows(focal_19nm, focal_21nm_TY_b, focal_22nm)
## 毛づくろい相手を表す列を追加
focal_nm_TY %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_nm_TY_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のTYとの毛づくろい時間を算出する。
focal_nm_TY_b %>%
filter(groom == "TY"|groom2 == "TY") %>%
group_by(subject) %>%
summarise(no_groom = n()) %>%
ungroup() -> no_groom_nm_TY最後に、毛づくろい時間割合を算出する。
3.2.4.2 近接時間割合
地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。
focal_nm_TY_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_nm_TY分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。
focal_prox_nm_TY %>%
group_by(subject) %>%
summarise(dur = n()) %>%
ungroup() -> prox_duration_nm_TYTYと3m以内に近接していたポイント数を算出する。
focal_prox_nm_TY %>%
filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>%
group_by(subject) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_nm_TY最後に、近接時間割合を算出する。
3.2.4.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_nm_TY <- mean(prop_groom_nm_TY$prop_groom)
## 平均近接時間割合
mean_prox_nm_TY <- mean(prop_prox_nm_TY$prop_prox)
## CSIの算出
prop_groom_nm_TY %>%
left_join(prop_prox_nm_TY, by = "subject") %>%
mutate(CSI_TY = ((prop_groom/mean_groom_nm_TY) + (prop_prox/mean_prox_nm_TY))/2) -> CSI_nm_TY算出した値が以下の通り。
CSI_nm_TY %>%
mutate(subject = fct_reorder(subject, CSI_TY)) %>%
ggplot(aes(x = subject, y = CSI_TY))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(CSI[i]))+
coord_flip(ylim = c(0,5.3)) +
scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
expand= c(0,0))+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 10.5),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(hjust = 0.5)) 
3.2.5 CSIの算出(調査期間ごと)
3.2.5.1 毛づくろい時間割合の算出
各メスの追跡時間(総瞬間サンプリングポイント数)を調査期間ごとに算出する。
focal_all_TY_b %>%
group_by(subject, study_period) %>%
summarise(dur = n()) %>%
ungroup() -> focal_duration_sp_TY各個体のTYとの毛づくろい時間を算出する。
focal_all_TY_b %>%
filter(groom == "TY"|groom2 == "TY") %>%
group_by(subject, study_period) %>%
summarise(no_groom = n()) %>%
ungroup() -> no_groom_sp_TY最後に、毛づくろい時間割合を算出する。
3.2.5.2 近接時間割合
地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。
調査期間ごとに分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。
focal_prox_TY %>%
group_by(subject, study_period) %>%
summarise(dur = n()) %>%
ungroup() -> prox_duration_sp_TYTYと3m以内に近接していたポイント数を算出する。
focal_prox_TY %>%
filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>%
group_by(subject, study_period) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_sp_TY最後に、近接時間割合を算出する。
3.2.5.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## CSIの算出
prop_groom_sp_TY %>%
left_join(prop_prox_sp_TY, by = c("subject","study_period")) %>%
group_by(study_period) %>%
mutate(mean_groom = mean(prop_groom),
mean_prox = mean(prop_prox)) %>%
ungroup() %>%
mutate(CSI_TY = ((prop_groom/mean_groom) + (prop_prox/mean_prox))/2) -> CSI_sp_TY算出した値が以下の通り。
CSI_sp_TY %>%
ggplot(aes(x = study_period, y = CSI_TY))+
geom_point()+
geom_line(aes(group = subject))+
theme_bw()+
facet_rep_wrap(~subject, repeat.tick.labels = TRUE)+
ylab(expression(CSI[i]))
3.3 ITとメスのCSIの算出
3.3.1 データの読み込みと加工
まずは、ITの出入りがあった2019年交尾期、2021年非交尾期についてTYが確認できた時間帯のデータを読み込む。
IT_19m <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
sheet = "male_presence_long") %>%
filter(maleID == "IT") %>%
mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
hour(.),minute(.))) %>%
mutate(date = as_date(date)) %>%
## 分派集団にいるときは除く
filter(groupID != 41 & groupID != 53 & groupID != 54)
IT_21nm <- read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
filter(maleID == "IT") %>%
mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
hour(.),minute(.))) %>%
mutate(date = as_date(date)) 上記2つの期間の個体追跡データに、各個体追跡セッション中にITが確認できていたかを表す列を追加し、ITが確認できたものだけを抽出する(focal_19m_b、focal_21nm_b)。
focal_19m %>%
distinct(no_focal, date, start_time, fin_time) %>%
left_join(IT_19m %>% select(date, male_presence, first, last), by = "date") %>%
## 個体追跡セッション中にITが確認されていたか
mutate(IT = ifelse((first <= start_time & last >= fin_time)|
(first >= start_time & first <= fin_time)|
(last >= start_time & last <= fin_time), 1, 0)) %>%
replace_na(list(IT = 0)) -> focal_19m_IT
focal_19m %>%
left_join(focal_19m_IT %>% select(no_focal, IT), by = "no_focal") %>%
filter(IT == "1") -> focal_19m_IT_b
focal_21nm %>%
distinct(no_focal, date, start_time, fin_time) %>%
left_join(IT_21nm %>% select(date, male_presence, first, last), by = "date") %>%
## 個体追跡セッション中にITが確認されていたか
mutate(IT = ifelse((first <= start_time & last >= fin_time)|
(first >= start_time & first <= fin_time)|
(last >= start_time & last <= fin_time), 1, 0)) %>%
replace_na(list(IT = 0)) -> focal_21nm_IT
focal_21nm %>%
left_join(focal_21nm_IT %>% select(no_focal, IT), by = "no_focal") %>%
filter(IT == "1") -> focal_21nm_IT_b交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_IT_b、focal_19m_IT_b)。
focal_18m %>%
left_join(female_18m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_18m_IT_b
focal_19m_IT_b %>%
left_join(female_19m %>%
rename(subject = femaleID) %>%
select(date, subject, rs2) , by = c("date", "subject")) %>%
filter(rs2 == "0") -> focal_19m_IT_b3.3.2 CSIの算出(全調査期間)
3.3.2.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 全期間を結合
focal_all_IT <- bind_rows(focal_18m_IT_b, focal_19m_IT_b, focal_19nm,
focal_21nm_IT_b)
## 毛づくろい相手を表す列を追加
focal_all_IT %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_IT_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のITとの毛づくろい時間を算出する。
focal_all_IT_b %>%
filter(groom == "IT"|groom2 == "IT") %>%
group_by(subject) %>%
summarise(no_groom = n())%>%
ungroup() -> no_groom_IT最後に、毛づくろい時間割合を算出する。
3.3.2.2 近接時間割合
地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。
focal_all_IT_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_IT分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。
ITと3m以内に近接していたポイント数を算出する。
focal_prox_IT %>%
filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>%
group_by(subject) %>%
summarise(no_prox = n())%>%
ungroup() -> no_prox_IT最後に、近接時間割合を算出する。
3.3.2.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_IT <- mean(prop_groom_IT$prop_groom)
## 平均近接時間割合
mean_prox_IT <- mean(prop_prox_IT$prop_prox)
## CSIの算出
prop_groom_IT %>%
left_join(prop_prox_IT, by = "subject") %>%
mutate(CSI_IT = ((prop_groom/mean_groom_IT) + (prop_prox/mean_prox_IT))/2) -> CSI_IT算出した値が以下の通り。Kitは一日ITに長時間毛づくろいした日があり、それが強く影響したようだ。
CSI_IT %>%
mutate(subject = fct_reorder(subject, CSI_IT)) %>%
ggplot(aes(x = subject, y = CSI_IT))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Times New Roman")+
xlab("")+
ylab(expression(paste(CSI, " with ", italic(IT))))+
labs(title = "B")+
coord_flip(ylim = c(0,6.0)) +
scale_y_continuous(breaks = seq(0,7, by= 0.5),
expand= c(0,0))+
theme_bw()+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "bold.italic",
hjust=0,
family = "Times New Roman"),
axis.text.x = element_text(family = "Times New Roman",
face = "bold"),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8,
family = "Times New Roman"),
axis.title.x = element_text(size = 10.5,
family = "Times New Roman"),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(family = "Times New Roman",
hjust = 0)) -> p_CSI_IT
p_CSI_IT
3.3.3 CSIの算出(2018交尾期～2019交尾期)
3.3.3.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 期間を結合
focal_fh_IT <- bind_rows(focal_18m_IT_b, focal_19m_IT_b, focal_19nm)
## 毛づくろい相手を表す列を追加
focal_fh_IT %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_IT_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のITとの毛づくろい時間を算出する。
focal_fh_IT_b %>%
filter(groom == "IT"|groom2 == "IT") %>%
group_by(subject) %>%
summarise(no_groom = n())%>%
ungroup() -> no_groom_fh_IT最後に、毛づくろい時間割合を算出する。
3.3.3.2 近接時間割合
地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。
focal_fh_IT_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_fh_IT分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。
ITと3m以内に近接していたポイント数を算出する。
focal_prox_fh_IT %>%
filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>%
group_by(subject) %>%
summarise(no_prox = n())%>%
ungroup() -> no_prox_fh_IT最後に、近接時間割合を算出する。
3.3.3.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_fh_IT <- mean(prop_groom_fh_IT$prop_groom)
## 平均近接時間割合
mean_prox_fh_IT <- mean(prop_prox_fh_IT$prop_prox)
## CSIの算出
prop_groom_fh_IT %>%
left_join(prop_prox_fh_IT, by = "subject") %>%
mutate(CSI_IT = ((prop_groom/mean_groom_fh_IT) + (prop_prox/mean_prox_fh_IT))/2) -> CSI_fh_IT算出した値が以下の通り。
CSI_fh_IT %>%
mutate(subject = fct_reorder(subject, CSI_IT)) %>%
ggplot(aes(x = subject, y = CSI_IT))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(CSI[i]))+
coord_flip(ylim = c(0,10)) +
scale_y_continuous(breaks = seq(0,10, by= 1),
expand= c(0,0))+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 10.5),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(hjust = 0.5)) 
3.3.4 CSIの算出(非交尾期のみ)
3.3.4.1 毛づくろい時間割合の算出
算出のためにデータを加工する。
## 全期間を結合
focal_nm_IT <- bind_rows(focal_19nm, focal_21nm_IT_b)
## 毛づくろい相手を表す列を追加
focal_nm_IT %>%
mutate(groom = ifelse(groomer == subject,groomee,groomer),
groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_nm_IT_b各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。
各個体のITとの毛づくろい時間を算出する。
focal_nm_IT_b %>%
filter(groom == "IT"|groom2 == "IT") %>%
group_by(subject) %>%
summarise(no_groom = n()) %>%
ungroup() -> no_groom_nm_IT最後に、毛づくろい時間割合を算出する。
3.3.4.2 近接時間割合
地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。
まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。
focal_nm_IT_b %>%
filter(activity %in% c("F","R","G") & T_G == "G") %>%
select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>%
replace_na(list(groom = "NA", groom2 = "NA")) %>%
filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_nm_IT分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。
focal_prox_nm_IT %>%
group_by(subject) %>%
summarise(dur = n()) %>%
ungroup() -> prox_duration_nm_ITITと3m以内に近接していたポイント数を算出する。
focal_prox_nm_IT %>%
filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>%
group_by(subject) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_nm_IT最後に、近接時間割合を算出する。
3.3.4.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。
## 平均毛づくろい時間割合
mean_groom_nm_IT <- mean(prop_groom_nm_IT$prop_groom)
## 平均近接時間割合
mean_prox_nm_IT <- mean(prop_prox_nm_IT$prop_prox)
## CSIの算出
prop_groom_nm_IT %>%
left_join(prop_prox_nm_IT, by = "subject") %>%
mutate(CSI_IT = ((prop_groom/mean_groom_nm_IT) + (prop_prox/mean_prox_nm_IT))/2) -> CSI_nm_IT算出した値が以下の通り。
CSI_nm_IT %>%
mutate(subject = fct_reorder(subject, CSI_IT)) %>%
ggplot(aes(x = subject, y = CSI_IT))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(CSI[i]))+
coord_flip(ylim = c(0,7)) +
scale_y_continuous(breaks = seq(0,10, by= 1),
expand= c(0,0))+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 10.5),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(hjust = 0.5)) 
3.3.5 CSIの算出(調査期間ごと)
3.3.5.1 毛づくろい時間割合の算出
各メスの追跡時間(総瞬間サンプリングポイント数)を調査期間ごとに算出する。
focal_all_IT_b %>%
group_by(subject, study_period) %>%
summarise(dur = n()) %>%
ungroup() -> focal_duration_sp_IT各個体のITとの毛づくろい時間を算出する。
focal_all_IT_b %>%
filter(groom == "IT"|groom2 == "IT") %>%
group_by(subject, study_period) %>%
summarise(no_groom = n()) %>%
ungroup() -> no_groom_sp_IT最後に、毛づくろい時間割合を算出する。
3.3.5.2 近接時間割合
地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。
調査期間ごとに分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。
focal_prox_IT %>%
group_by(subject, study_period) %>%
summarise(dur = n()) %>%
ungroup() -> prox_duration_sp_ITITと3m以内に近接していたポイント数を算出する。
focal_prox_IT %>%
filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>%
group_by(subject, study_period) %>%
summarise(no_prox = n()) %>%
ungroup() -> no_prox_sp_IT最後に、近接時間割合を算出する。
3.3.5.3 CSIの算出
それでは、上記のデータを用いてCSIを算出する。なお、平均毛づくろい時間割合が0の調査期間は分子の左側を1とした。
## CSIの算出
prop_groom_sp_IT %>%
left_join(prop_prox_sp_IT, by = c("subject","study_period")) %>%
group_by(study_period) %>%
mutate(mean_groom = mean(prop_groom),
mean_prox = mean(prop_prox)) %>%
ungroup() %>%
mutate(CSI_IT = ifelse(mean_groom != 0,
((prop_groom/mean_groom) + (prop_prox/mean_prox))/2,
(1 + (prop_prox/mean_prox))/2)) -> CSI_sp_IT算出した値が以下の通り。
CSI_sp_IT %>%
ggplot(aes(x = study_period, y = CSI_IT))+
geom_point()+
geom_line(aes(group = subject))+
theme_bw()+
facet_rep_wrap(~subject, repeat.tick.labels = TRUE)+
ylab(expression(CSI[i]))
3.4 分母をTYとIT両方の頻度の平均値にする場合
CSIを算出する際の分母を、TYとITのデータの両方を合わせた平均にします。
3.4.1 CSIの算出
算出は以下のように行えます。
## 平均毛づくろい時間割合
mean_groom_TYIT <- mean(c(prop_groom_TY$prop_groom,
prop_groom_IT$prop_groom))
## 平均近接時間割合
mean_prox_TYIT <- mean(c(prop_prox_TY$prop_prox,
prop_prox_IT$prop_prox))
## CSIの算出
### 平均値
prop_groom_TY %>%
left_join(prop_prox_TY, by = "subject") %>%
mutate(CSI_TY = ((prop_groom/mean_groom_TYIT) + (prop_prox/mean_prox_TYIT))/2) -> CSI_TY_combined
prop_groom_IT %>%
left_join(prop_prox_IT, by = "subject") %>%
mutate(CSI_IT = ((prop_groom/mean_groom_TYIT) + (prop_prox/mean_prox_TYIT))/2) -> CSI_IT_combined 3.4.2 図示
図示すると以下のようになります。
CSI_TY_combined %>%
mutate(subject = fct_reorder(subject, CSI_TY)) %>%
ggplot(aes(x = subject, y = CSI_TY))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(paste(CSI, " with ", italic(TY))))+
labs(title = "(a)")+
coord_flip(ylim = c(0,5.15)) +
scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
expand= c(0,0))+
theme_bw()+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0,
family = "Times New Roman"),
axis.text.x = element_text(family = "Times New Roman"),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8,
family = "Times New Roman"),
axis.title.x = element_text(size = 10.5,
family = "Times New Roman"),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(family = "Times New Roman",
hjust = 0)) -> p_CSI_TY_combined
# ggsave("figures/p_CSI_TY_combined.tif",
# p_CSI_TY_combined, width = 100, height = 100, dpi = 2000,
# units = "mm")
CSI_IT_combined %>%
mutate(subject = fct_reorder(subject, CSI_IT)) %>%
ggplot(aes(x = subject, y = CSI_IT))+
geom_col(color = "black", fill = "white", size =0.3)+
theme_classic(base_size=12, base_family = "Arial")+
xlab("")+
ylab(expression(paste(CSI, " with ", italic(IT))))+
labs(title = "(b)")+
coord_flip(ylim = c(0,4)) +
scale_y_continuous(breaks = seq(0,4, by= 0.5),
expand= c(0,0))+
theme_bw()+
theme(aspect.ratio=1,
axis.text.y = element_text(face = "italic",
hjust=0,
family = "Times New Roman"),
axis.text.x = element_text(family = "Times New Roman"),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 8,
family = "Times New Roman"),
axis.title.x = element_text(size = 10.5,
family = "Times New Roman"),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
plot.title = element_text(family = "Times New Roman",
hjust = 0)) -> p_CSI_IT_combined
# ggsave("figures/p_CSI_IT_combined.tif",
# p_CSI_IT_combined, width = 100, height = 100, dpi = 2000,
# units = "mm")
p_CSI_TYIT_combined <- p_CSI_TY_combined + p_CSI_IT_combined
# ggsave("figures/p_CSI_TYIT_combined.tif", p_CSI_TYIT_combined,
# width = 180, height = 105, units = "mm", dpi = 2000)3.5 算出に使用したメスの個体追跡データの概要
以下に、全期間のCSIを算出するのに使用した各メスの個体追跡データの総時間(総瞬間サンプリングポイント数)の情報をまとめる。
focal_duration_TY %>%
left_join(focal_duration_IT, by = "subject") %>%
rename(femaleID = 1, "Total points(TY)" = 2, "Total points (IT)" = 3) %>%
flextable() %>%
colformat_double(digits=2) %>%
set_table_properties(layout="autofit",width = 1) %>%
autofit(add_w = 0.2) %>%
font(part = "header", fontname = "Times New Roman") %>%
font(part = "body", j=1:3, fontname = "Times New Roman") %>%
theme_zebra() %>%
hline_bottom() %>%
hline_top() %>%
align(j=2:3,part = "all",align = "center") -> table_focal_duration
table_focal_durationfemaleID | Total points(TY) | Total points (IT) |
|---|---|---|
Aka | 2,306 | 1,575 |
Ako | 2,459 | 1,704 |
Hen | 2,621 | 1,689 |
Hot | 2,268 | 1,603 |
Kil | 2,561 | 1,655 |
Kit | 2,302 | 1,545 |
Koh | 2,401 | 1,731 |
Kol | 2,627 | 1,715 |
Kun | 2,099 | 1,344 |
Kur | 1,120 | 1,120 |
Mal | 2,458 | 1,640 |
Mei | 2,742 | 1,727 |
Mif | 1,225 | 547 |
Mik | 2,413 | 1,717 |
Ntr | 2,800 | 1,801 |
Tam | 1,458 | 1,458 |
Ten | 2,618 | 1,726 |
Tot | 2,720 | 1,848 |
4 個体追跡データの加工
以下では、個体追跡データを分析のために加工します。
4.2 データの加工
4.2.1 メスの情報
まず、個体追跡データに以下のような加工を行い、新たに列を作成します。
groupID: 追跡メスがどの集団にいたか
max_female: 最大メス数
no_female: 確認メス数
no_est: 発情メス数
## 各フォーカル個体がどの集団で観察されたか
female_presence_group <- group_all %>%
pivot_longer(Kur:Yun,
names_to = "subject",
values_to = "presence") %>%
filter(presence == 1) %>%
## 2020年11月13の2つ目の集団は除く
filter(!groupID %in% c("m20_52")) %>%
select(date, groupID, subject)
## 各観察日の最大個体
max_female <- female_time %>%
pivot_longer(cols = Kur:Cur,
names_to = "femaleID",
values_to = "presence") %>%
filter(is.na(presence)|(presence != "DD" & presence != "NS")) %>%
group_by(date) %>%
summarise(max_female = n())
## 各観察日の確認メス数(6歳以上)と発情メス数
no_female <- group_all %>%
select(groupID, study_period, date, Kur:Yun) %>%
select(-c(TY,IT, LK, KR, KM, TG)) %>%
pivot_longer(cols = Kur:Yun,
names_to = "femaleID",
values_to = "presence") %>%
left_join(att) %>%
## 6歳以上の個体のみを抽出
filter(age >= 6) %>%
left_join(female_all %>% select(date, femaleID, rs2)) %>%
## groupIDごとに個体数を算出
group_by(date, groupID, study_period) %>%
summarise(no_female = sum(presence, na.rm = TRUE),
no_est = sum(rs2, na.rm = TRUE)) %>%
ungroup()
## 結合
focal_combined_b <- focal_combined %>%
mutate(date = as_date(date)) %>%
left_join(female_presence_group,
by = c("subject", "date")) %>%
left_join(max_female,
by = c("date")) %>%
left_join(no_female)4.2.2 オスの情報
続いて、交尾期については追跡日に確認された群れ外オス数、オス数を追加します。
まず、生データを用いて群れ外オス数とオス数を算出します。
no_ntm_m18 <- male_18m %>%
mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK"),
0, 1)) %>%
group_by(date) %>%
summarise(no_male = sum(presence, na.rm = TRUE),
no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))
no_ntm_m19 <- male_19m %>%
mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK"),
0, 1)) %>%
group_by(date) %>%
summarise(no_male = sum(presence, na.rm = TRUE),
no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))
no_ntm_m20 <- male_20m %>%
mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK", "KM"),
0, 1)) %>%
group_by(date) %>%
summarise(no_male = sum(presence, na.rm = TRUE),
no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))
no_ntm_m21 <- male_21m %>%
mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK", "KM", "TG"),
0, 1)) %>%
group_by(date) %>%
summarise(no_male = sum(presence, na.rm = TRUE),
no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))
sum_ntm <- bind_rows(no_ntm_m18, no_ntm_m19, no_ntm_m20, no_ntm_m21)これを結合します。
4.2.3 TYとITの確認状況
調査期間中、第一位オスのタイヨウ(TY)と第二位オスのイツモ(IT)の群れへの出入りが頻繁に観察されました。
そこで、彼らが個体追跡時にいたか否かについての列(TY、IT)を追加します。
まず、彼らが確認できた時刻に関するデータを読み込む。
TYIT_presence_time <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "m19") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last) %>%
bind_rows(read_excel("../Data/data/2020mating/2020mating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "m20") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
bind_rows(read_excel("../Data/data/2021mating/2021mating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "m21") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
bind_rows(read_excel("../Data/data/2019nonmating/2019nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "nm19") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
bind_rows(read_excel("../Data/data/2020nonmating/2020nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "nm20") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
bind_rows(read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "nm21") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
bind_rows(read_excel("../Data/data/2022nonmating/2022nonmating_raw.xlsx",
sheet = "male_presence_long") %>%
mutate(study_period = "nm22") %>%
mutate(groupID = str_c(study_period,"_", groupID)) %>%
select(maleID, date, groupID, male_presence, first, last)) %>%
mutate(across(c(first, last),
~make_datetime(year(date), month(date), mday(date), hour(.),minute(.))))
## 横長にする
TYIT_presence_time_wide <- TYIT_presence_time %>%
filter(maleID == "TY") %>%
rename(first_TY = first,
last_TY = last,
TY = male_presence) %>%
select(-maleID) %>%
left_join(TYIT_presence_time %>%
filter(maleID == "IT") %>%
rename(first_IT = first,
last_IT = last,
IT = male_presence) %>%
select(-maleID)) %>%
mutate(date = as_date(date))
TYIT_presence_time_wideこれを利用し、個体追跡中にTYまたはITがいたと考えられるかを表す列を作成します。
TY_presence/IT_presence1: 個体追跡開始前または個体追跡中にTYを確認、かつTYがいなくなったと思われる前に開始
0: その他
focal_combined_c %>%
left_join(group_all %>% select(date, groupID, TY, IT)) %>%
left_join(TYIT_presence_time_wide %>% select(-TY, -IT)) %>%
mutate(TY_presence = case_when(
study_period == "m18" ~ TY,
is.na(first_TY) ~ 0,
start_time >= first_TY & start_time <= last_TY ~ 1,
fin_time >= first_TY ~ 1,
start_time >= last_TY ~ 0,
fin_time <= first_TY ~ 0,
.default = 1
)) %>%
mutate(IT_presence = case_when(
study_period == "m18" ~ IT,
is.na(first_IT) ~ 0,
start_time >= first_IT & start_time <= last_IT ~ 1,
fin_time >= first_IT ~ 1,
start_time >= last_IT ~ 0,
fin_time <= first_IT ~ 0,
.default = 1
)) -> focal_combined_d4.2.4 フォーカルを攻撃したオスのID
最後に、フォーカルが攻撃を受けた際の攻撃したオスのIDを表す列を作成します。
focal_combined_d %>%
## NAになっているところを0に
replace_na(list(agg_focal = 0)) %>%
replace_na(list(victim1 = "NA", victim2 = "NA", victim3 = "NA",
aggressor1 = "NA", aggressor2 = "NA", aggressor = "NA")) %>%
mutate(if_victim1 = ifelse(agg_focal == 1 & str_detect(victim1, subject), 1, 0),
if_victim2 = ifelse(agg_focal == 1 & str_detect(victim2, subject), 1, 0),
if_victim3 = ifelse(agg_focal == 1 & str_detect(victim3, subject), 1, 0)) %>%
mutate(aggressor_focal = case_when(
if_victim1 == 1 & if_victim2 == 1 & if_victim3 == 1 ~ str_c(aggressor1,aggressor2,aggressor3,sep =","),
if_victim1 == 1 & if_victim2 == 1 ~ str_c(aggressor1, aggressor2,sep =","),
if_victim1 == 1 & if_victim3 == 1 ~ str_c(aggressor1, aggressor3,sep =","),
if_victim2 == 1 & if_victim3 == 1 ~ str_c(aggressor2, aggressor3,sep =","),
if_victim1 == 1 ~ aggressor1,
if_victim2 == 1 ~ aggressor2,
if_victim3 == 1 ~ aggressor3,
.default = NA
)) %>%
## aggressor2には、追跡個体が同時刻に2頭から攻撃を受けたときのみ値が入る
separate(aggressor_focal, into = c("aggressor_focal1","aggressor_focal2"),
sep = ",") -> focal_combined_final加工したデータは以下の通り。